MODULE fefmm_quadrature

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_get_roots_and_weights

CONTAINS

!-------------------------------------------------------------------------------
! Gauss-Legendre 10-point quadrature over hierarchical cells

   SUBROUTINE fefmm_get_roots_and_weights_old(npoints,tmax,roots,weights)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: npoints
      REAL(REALK),   INTENT(IN)  :: tmax
      REAL(REALK),   INTENT(OUT) :: roots(:), weights(:)

      REAL(REALK), PARAMETER :: x(5) = (/0.1488743389d0, 0.4333953941d0, &
                                         0.6794095682d0, 0.8650633666d0, &
                                         0.9739065285d0/)
      REAL(REALK), PARAMETER :: wt(5) = (/0.2955242247d0, 0.2692667193d0, &
                                          0.2190863625d0, 0.1494513491d0, &
                                          0.0666713443d0/)
      REAL(REALK)   :: t1,t2
      REAL(REALK)   :: xr,xm,dx
      INTEGER(INTK) :: tn, np,nc, ncells, k

      ncells = INT(npoints/10)

      k = 0
      DO nc = 1, ncells

         t2 = tmax*(0.5d0**(nc-1))
         t1 = t2*0.5d0

         IF (nc == ncells) t1 = 0d0
         xm = 0.5d0*(t2+t1)
         xr = 0.5d0*(t2-t1)

         DO np = 1, 5
            k = k+2
            dx = xr*x(np)

            roots(k)     = xm+dx
            roots(k-1)   = xm-dx
            weights(k)   = wt(np)*xr
            weights(k-1) = wt(np)*xr

         END DO
      END DO

   END SUBROUTINE fefmm_get_roots_and_weights_old

!-------------------------------------------------------------------------------
! General routine to return roots and weights of the Gauss-Legendre
! n-point quadrature formula, given lower,x1 and upper,x2 integration limits

   SUBROUTINE fefmm_get_roots_and_weights(n,x2,roots,weights)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: n
      REAL(REALK),   INTENT(IN)  :: x2
      REAL(REALK),   INTENT(OUT) :: roots(:), weights(:)

      REAL(REALK), PARAMETER :: x1 = zero
      REAL(REALK)   :: xm,xl, z,z1, pp,p3,p2,p1
      INTEGER(INTK) :: m,i,j

      m = (n+1)/2
      xm = half*(x2+x1)
      xl = half*(x2-x1)

      DO i = 1, m
         z = COS(PI*(i-0.25d0)/(n+half))

         DO ! until z converged
            p1 = one
            p2 = zero
            DO j = 1, n
               p3 = p2
               p2 = p1
               p1 = ((two*j-one)*z*p2 - (j-one)*p3) / j
            END DO
            pp = n*(z*p1-p2) / (z*z-one) 
            z1 = z
            z = z1 - p1/pp
            IF ( ABS(z-z1) < 1d-15) EXIT
         END DO

         roots(i) = xm - xl*z
         roots(n+1-i) = xm + xl*z
         weights(i) = two*xl / ((one-z*z)*pp*pp)
         weights(n+1-i) = weights(i)

      END DO

   END SUBROUTINE fefmm_get_roots_and_weights

!-------------------------------------------------------------------------------

END MODULE fefmm_quadrature
